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Abstract 

The well-scaled transition to the diffusion limit in the framework 
of the theory of continuous-time random walk (CTRW) is presented 
starting from its representation as an infinite series that points out the 
subordinated character of the CTRW itself. We treat the CTRW as 
a combination of a random walk on the axis of physical time with a 
random walk in space, both walks happening in discrete operational 
time. In the continuum limit we obtain a (generally non-Markovian) 
diffusion process governed by a space-time fractional diffusion equa- 
tion. The essential assumption is that the probabilities for waiting 
times and jump-widths behave asymptotically like powers with nega- 
tive exponents related to the orders of the fractional derivatives. By 
what we call parametric subordination, applied to a combination of a 
Markov process with a positively oriented Levy process, we generate 
and display sample paths for some special cases. 
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1 Introduction 



Surveying the literature of the past 15 years we can observe an ever increas- 
ing interest in modelling anomalous diffusion processes, namely in diffusion 
processes deviating essentially from Gaussian behaviour which is charac- 
terized by evolution of the second centered moment like the first power of 
time. The reader interested to these processes is referred to several educa- 
tional/review papers and books, including [21 HZl HSl 113 I2H EH EH E3 

S2i sal esi esi eh eu 

In Section 2, we recall the simplest models for anomalous diffusion based 
on fractional calculus. They are obtained by replacing in the classical dif- 
fusion equation the partial derivatives with respect to space and/or time 
by derivatives of non- integer order, in such a way that the resulting Green 
function can still be interpreted as a probability density evolving in time 
differently from the Gaussian type. 

A more general approach to anomalous diffusion is provided by the so- 
called continuous time random walk (CTRW) introduced in Statistical Me- 
chanics by Montroll and Weiss [39] , see also [371 EH US [60] , which differs 
from the usual models in that the steps of the walker occur at random times 
generated by a renewal process. The sojourn probability density of this pro- 
cess is known to be governed by an integral equation and expressed in terms 
of a relevant series expansion, as it will be recalled in Section 3. The concept 
of CTRW, can be understood by considering a random walk subordinated 
to a renewal process, see e.g. [9], as pointed out by a number of authors, see 

e.g. [n m Eg Eaisa Elic- 
it is well known that the space-time fractional diffusion (STFD) equation 

and its variants, including the fractional Fokker-Planck equation, can be 

derived from the CTRW integral equation, see e.g. [H El EQl [Ml USE [361 

021 EH E31 El], and references therein. More rigorously the passage from 

CTRW to STFD can be carried out via a properly scaled transition to the 

diffusion limit (under appropriate assumptions on waiting times and jumps), 

as shown in [58] and in a number of papers of our research group, see e.g. 

[m mug gg. 

In this paper, we offer another scheme of well-scaled transition to the 
diffusion limit, a scheme based on a modified concept of subordination that 
we call parametric subordination. To this purpose we lay open our general 
view of subordination in stochastic processes in Section 4. Then in Sec- 
tion 5, starting from the series expansion of the sojourn probability density 
in CTRW we arrive in a well-scaled limit process at the relevant integral 
formula for subordination. 
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Finally, we consider the problem of how to construct the sample paths 
for the STFD based on the above diffusion limit of the CTRW. In Section 6, 
we explain what we mean by parametric subordination whereas in Section 
7, we describe the numerical procedure and provide sample paths for four 
case studies. The main conclusions are drawn in Section 8. 

2 The space-time fractional diffusion 

We begin by considering the Cauchy problem for the (spatially one- 
dimensional) space-time fractional diffusion equation 

t D^u(x,t) = x D?u(x,t), u(x,0) = S(x), i£R, t>0, (2.1) 

where {a , , f3} are real parameters restricted to the ranges 

0<a<2, \6\ <mm{a,2 - a} , < (3 < 1 . (2.2) 

Here tD* denotes the Caputo fractional derivative of order f3, acting on the 
time variable t, and x Dg denotes the Riesz-Feller fractional derivative of 
order a and skewness 9, acting on the space variable x. Let us note that the 
solution u(x,t) of the Cauchy problem (2.1), known as the Green function 
or fundamental solution of the space-time fractional diffusion equation, is a 
probability density in the spatial variable x, evolving in time t. In the case 
a = 2 and (5 = 1 we recover the standard diffusion equation for which the 
fundamental solution is the Gaussian density with variance a 2 = 2t. 

Writing, with Re[s] > do, k € R, the transforms of Laplace and Fourier 

as 

C{f(t);s} = f(s):= / e- st f(t)dt, 
Jo 

/+oo . 
e lKX g{x) dx, 
-oo 

we have the corresponding transforms of tD* f(t) and xD^g^x) as 

£ { tD? f(t)} = sf 3 fis)- S f 3 - 1 fiO) , (2.3) 

F{ x Dfg(x)} = -\K\ a i es[ % nK gi K ) . (2.4) 

Notice that ^sig nK = exp[z (sign k) 6 vr/2]. For the mathematical details 
the interested reader is referred to [131 113 HI] on t ne Caputo derivative, and 
to [44j on the Feller potentials. For the general theory of pseudo-differential 
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operators and related Markov processes the interested reader is referred to 
the excellent volumes by Jacob |22j . 

For our purposes let us here confine ourselves to recall the representation 
in the Laplace- Fourier domain of the (fundamental) solution of (2.1) as it 
results from the application of the transforms of Laplace and Fourier. Using 
S(k) = 1 we have from (2.1) 

s $ 5(«, s)- s f 3 - 1 = -\ K \ a i e si e n K a) , 

hence 

u(k,s) = —q . (2.5) 

s p + \ K \ assign K 

For explicit expressions and plots of the fundamental solution of (2.1) in the 
space-time domain we refer the reader to [27] . There, starting from the fact 
that the Fourier transform u(k, t) can be written as a Mittag-Lefner function 
with complex argument, the authors have derived a Mellin-Barnes integral 
representation of u(x, t) with which they have proved the non-negativity 
of the solution for values of the parameters {a, 9, /?} in the range (2.2) 
and analyzed the evolution in time of its moments. In particular for {0 < 
a<2,/3 = l}we obtain the stable densities of order a and skewness 9. 
The representation of u(x, t) in terms of Fox //-functions can be found in 
[29] . We note, however, that the solution of the STFD Equation (2.1) and 
its variants has been investigated by several authors as pointed out in the 
bibliography in [27] : here we refer to some of them, [H El l33l [35] , where the 
connection with the CTRW was also pointed out. 



3 The continuous-time random walk 

The name continuous time random walk (CTRW) became popular in physics 
after Montroll, Weiss and Scher (just to cite the pioneers) in the 1960s and 
1970s published a celebrated series of papers on random walks for modelling 
diffusion processes on lattices, see e.g. [37J, [39] , and the book by Weiss [60] 
with references therein. CTRWs are rather good and general phenomenolog- 
ical models for diffusion, including processes of anomalous transport, that 
can be understood in the framework of the classical renewal theory, as stated 
e.g. in the booklet by Cox [9]. In fact a CTRW can be considered as a com- 
pound renewal process (a simple renewal process with reward) or a random 
walk subordinated to a simple renewal process. 

Basic notions of the CTRW theory, that hereafter we briefly recall for 
the readers' convenience, are the master equation (in integral form) for the 
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sojourn probability density, its Fourier-Laplace representation (known as 
the Montroll- Weiss formula) and its series representation. 

A CTRW is generated by a sequence of independent identically dis- 
tributed (iid) positive random waiting times T\, T2, T3, . . . , each having the 
same probability density function <j)(t) , t > , and a sequence of iid random 
jumps X\ , X2,X%, . . . , in R , each having the same probability density w(x) , 
x £ R. 

Let us remark that, for ease of language, we use the word density also 
for generalized functions in the sense of Gel'fand and Shilov [llj, that can 
be interpreted as probability measures. Usually the probability density func- 
tions are abbreviated by pdf. We recall that (j)(t) > with f Q °° <f)(t) dt = 1 
and w(x) > with w(x) dx = 1. 

Setting to = , t n = T\ + T 2 + . . . T n for n G N , the wandering particle 
makes a jump of length X n in instant t n , so that its position is xq = 
for < t < T\ = t\ , and x n = X\ + X2 + . . . X n , for t n < t < t n+ \ . 
We require the distribution of the waiting times and that of the jumps to 
be independent of each other. So, we have a compound renewal process (a 
renewal process with reward), compare [9j. 

By natural probabilistic arguments we arrive at the integral equation for 
the probability density p(x, t) (a density with respect to the variable x) of 
the particle being in point x at instant t , see e.g. [HI [161 EHl H3 H71 135] . 



p(z,t) = 6(z)V(t) + f 4>{t-t') 
Jo 



w(x — x') p(x', t') dx' 



dt' , (3.1) 



in which the survival function 

q>(t) = I (P(t')dt' (3.2) 



denotes the probability that at instant t the particle is still sitting in its 
starting position x = . Clearly, (3.1) satisfies the initial condition p(x, 0) = 
5(x). In the Laplace- Fourier domain Eq. (3.1) reads as 

p(K, s) = ifr(s) + w(k) 4>(s)p(k,, s) , 

and using ^(s) = (1 — 0(s))/s , explicitly 

Kn,s)= 1 ~^ J„ . (3.3) 
s 1 — w(k) <j)(s) 

This Laplace-Fourier representation is known in physics as the the Montroll- 
Weiss equation, so named after the authors, see [3U], who derive it in 1965 
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as the basic equation for the CTRW. By inverting the transforms one can 
find the evolution p(x, t) of the sojourn density for time t running from zero 
to infinity. In fact, recalling that < 1 and \4>(s)\ < 1, if k 7^ and 

s / 0, Eq. (3.3) becomes 

00 00 
f( K ,s) = *(s) ^(s)w(K)] n = Y,Ms)w n (K) , (3.4) 

n=0 n=0 

and we promptly obtain the series representation of the continuous time 
random walk, see e.g. [9J (Ch. 8, Eq. (4)) or [60] (Eq.(2.101)), 

00 00 

p(x, t) = J2 V n(t) w n(x) = ®(t) 5(x) + ^ Mt) ™n(x) , (3-5) 
n=0 n=l 

where the functions v n (t) and w n (x) are obtained by repeated convolutions 
in time and in space, v n (t) = (V& * (p* n )(t), and w n (x) = (w* n )(x), re- 
spectively. In particular, v (t) = (<£ * 5)(t) = ^(t), vi(t) = (^ * 4>)(t), 
wq(x) = 5(x), w\{x) = w{x). In the R.H.S of Eq (3.5) we have isolated the 
first singular term related to the initial condition p{x, 0) = ^(0) 5(x) = 5(x). 
The representation (3.5) can be found without detour over (3.1) by direct 
probabilistic reasoning and transparently exhibits the CTRW as a subordi- 
nation of a random walk to a renewal process: it can be used as starting 
point to derive the Montroll- Weiss equation, as it was originally recognized 
by Montroll and Weiss [39]. Though (3.5), while being an attractive general 
formula, is unlikely to lead to explicit answers to rather simple problems, 
we consider it as a basic and useful formula for our analysis, as it will be 
shown later on. 

A special case of the integral equation (3.1) is obtained for the compound 
Poisson process where (p(t) = me~ mt (with some positive constant m). Then, 
the corresponding master equation reduces after some manipulations, that 
best are carried out in the Laplace-Fourier domain, to the Kolmogorov- Feller 
equation: 

— p(x, t) = — mp(x, t) + m / w(x — x') p(x' , t) dx . (3.6) 

Ot J —oo 

Then, the solution obtained via the series representation reads 

p(x,t)=Y,^-^ mt Mx). (3.7) 

k=0 

Note that only in this case the corresponding stochastic process is Marko- 
vian. 
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4 Subordination in stochastic processes 



In recent years a number of papers have appeared where explicitly or im- 
plicitly subordinated stochastic processes have been treated in view of their 
relevance in physical and financial applications, see e.g. [D El [32l EH1 
HHJ [501 EH [52l E3 ESI E2] and references therein. Historically, the notion of 
subordination was originated by Bochner, see 0[8]. 

We obtain the process X(t) of our proper interest in the form X(t) = 
Y(T*(t)) by randomizing the time clock of a stochastic process Y{t*) using a 
new clock t = T(t*), the non-decreasing right-continuous random functions 
t = T(t*) and t* = T*(t) being inverse (in the appropriate sense) to each 
other. The resulting process x = X{t) is said to be subordinated to the so- 
called parent process Y(t*), and i* is commonly referred to as the operational 
time. 

Our essential process for randomizing time is the process t = 
called by us the leading process. Our view is in contrast to that of Bochner's 
subordination adopted by Feller [10] and others, see e.g. [32J [57], who put 
into the foreground the inverse process i* = T*(t), which actually is a hit- 
ting time or first passage process, and after Feller often called the directing 
process. 

In particular, assuming Y(£*) to be a Markov process with a spa- 
tial probability density function (pdf) of x, evolving in operation time £*, 
qt. t {x) = q(x,t*), and T*(i) to be a process with non-negative, not neces- 
sarily independent, increments with pdf of i* depending on a parameter t, 
rt(t*) = r(£*,i), then the subordinated process X(t) = Y(T*(t)) is governed 
by the spatial pdf of x evolving with t, pt(x) = p(x, t), given by the integral 
formula of subordination (compare with Eq (7.1), Ch. X in [10J and with 



J o 

If the parent process Y(t*) is self-similar of the kind that its pdf qt,{x) is 
such that, with a probability density q{x) and a positive number 7, 



Eq. (3.1) in [28J) 




(4.1) 




(4.2) 



then Eq. (4.1) reads 




(4.3) 
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Remark : Feller [TO] and several other mathematicians, e.g. [22\ 145]. in their 
treatment of subordination, are mainly interested in Markov processes. Af- 
ter explicitly saying (in his Section X.7) that the subordinated process may 
happen to be non-Markovian, Feller immediately turns his attention to the 
search for conditions to be imposed on the directing process that ensure 
the subordinated process to be Markovian like the parent process. For our 
processes (see next Section) these conditions are in general not fulfilled. 

5 Subordination in continuous time random walk 

In fractional diffusion an intuitive understanding can be gained by formaliz- 
ing the transition from the series representations (3.4) and (3.5) of a general 
continuous time random walk (CTRW), known in the mathematical litera- 
ture as a renewal process with reward. We cannot survey the rich literature 
on the subject, but let us call here the reader's special attention to the most 
recent papers by Piryatinska, Saichev and Woyczynski [12] and Sokolov and 
Klafter [54] . These authors show in differing ways how fractional diffusion 
can be obtained from continuous time random walk. In contrast to these au- 
thors we lay out in all details our method of well-scaled transition to the dif- 
fusion limit, making explicit the meaning of long-time wide-space behaviour. 
For the general principle of well-scaledness we refer to [12} [T4"l \TE[ 08] . 

In the series representation (3.5) for the CTRW the running index n cor- 
responds to the so-called "operational time" t* in the subordination formula 
for a continuous (stable) process. We will pass in (3.5) to the diffusion limit 
under the "power law" assumptions (in the Laplace- Fourier domain) 

l-0(s)~As /3 , A>0, s^0+, (5.1) 

1-w(k) ~Ml«r^ SignK 5 M>0, k^O, (5.2) 

where /?, a and 6 are restricted as in (2.2). If < /3 < 1 and < a < 2, Eqs. 
(5.1) and (5.2) imply fat (power-law) tails for the densities 4>(t) and w(x); 
otherwise, for (3 = 1, Eq. (5.1) implies that (f>(t) has a finite first moment 
(e.g. the exponential pdf), and, for a = 2, Eq. (5.2) implies that w(x) has 
a finite second moment (e.g. the Gaussian pdf). For details we refer e.g. to 

The idea is to treat the series expansion (starting from n = 0) in (3.5) 
as an approximation to an improper Riemann integral. Being interested on 
behaviour in large time and wide space we change the units of measurement 
in order to make large time intervals and space distances appear numerically 
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of moderate size, moderate time intervals and space distances of small size. 
To this aim we replace waiting times T by tT, jumps X by hX, and then 
send the positive scaling factors r and h to zero, observing a scaling relation 
that will become mandatory in our calculations. For conciseness of our 
presentation we skip the analytical subtleties of interchanges of summations 
and integrations. For a strictly analytical derivation of our final integral 
equation of subordination we recommend [28]. 

For the CTRW this means replacing <j)(t) by <fi T (t) = (f>{t/r)/r, w(x) 
by u>h(x) = w(x/h)/h, correspondingly cp(s) by 4> T (s) = 4>(rs), w(k) by 
Wh( K ) = w(hn). Decorating (3.5) by indices h and r gives 



oo 



Ph,r(x, t)=} j V T ,n{t) Wh,n{x) , (5.3) 
n=0 

yielding in the Fourier-Laplace domain 

% >T (n,s) = f; 1 ~ ^ TS) U(rs)) n (w(hK)T • (5.4) 



n=0 

Separately we treat the powers (cj)(Ts)j and (w(/iK)) n , so avoiding the 
problematic simultaneous inversion of the diffusion limit from the Fourier- 
Laplace domain into the physical domain. 
Observing from (5.1) 

' (r S )) n ~(l-A(r S f) n , (5.5) 

we relate the running index n to the presumed operational time i* by 

n~^, (5.6) 

and for fixed s (as required by the continuity theorem of probability theory) , 
by sending r — > we get 

(raj) n ~ f 1 - A tV) U/{XT ' 3) exp (-U a?) . (5.7) 



Here s corresponds to physical time t, and in Laplace inversion we must 
treat i* as a parameter. Hence, in physical time exp(— i*s^) corresponds to 

gp(t,U) = t- 1/ Pgp{t- l/f3 t), (5.8) 
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with gp(s) = exp(— s 13 ). Here gp(t,t*) is the totally positively skewed sta- 
ble density (with respect to the variable t) evolving in operational time t* 
according to the "space"- fractional equation 

d 

— gp(t, U) = t D% gp(t, U) , g p (t, 0) = 6{t) , (5.9) 

where t is playing the role of the spatial variable. Analogously, observing 
from (5.2) 

{w{hn)) n ~ (l- fi(h\K\) a i esi Z nK y , (5.10) 
and with the aim of obtaining a meaningful limit we now set 

lih a ' y ' 

and find, by sending h — > + , the relation 

{w(hK)) n ~ (l - M%ir^ 8ignK ) W( ' l/ia) - exp (-^af^"^) , 

(5.12) 

the Fourier transform of a ^-skewed a-stable density f a ,e(x,t*) evolving in 
operational time t*. This density is the solution of the space-fractional 
equation 

d 

^-/a,0OM*) = xDg f a ,e{x,Q, f a fi{x, 0) = S(x) . (5.13) 

The two relations (5.6) and (5.11) between the running index n and the 
presumed operational time t* require the (asymptotic) scaling relation 

Xr^fih , (5.14) 

that for purpose of computation we simplify to 

\ T P = fih a . (5.15) 

Replacing t* by i* jn = nXr 13 , using the asymptotic results (5.7) and 
(5.12) obtained for the powers (0(rs)) and (w(hK)) n , furthermore noting 



1 ~ <P(ts) fl 
~ A T , 
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we finally obtain from (5.4) the Riemann sum (with increment At' 3 ) 

oo 
n=0 

and hence the integral 

Ph,A K i s ) ~ 5/3-1 / ex P 
j 

For the limiting process up(x,t) this means 

up(K,s) = / s^ _1 exp 

Observe that the RHS of this equation is just another way of writing the 
RHS of equation (2.5) which is the Laplace- Fourier solution of the STFD 
equation (2.1). By inverting the transforms we get after some manipulations 
(compare [32]) in physical space-time the integral formula of subordination 

/•oo 

up(x,t) = f a ,e(x,U) gp(U,t)dt* (5.19) 
J 

with 

gp(t*,t) = ^g/3 u l/!3 - 1 (5.20) 

standing for the density r t (t*) in equation (4.1). 

There are two processes involved. One is the unidirectional motion along 
the i* axis representing the operational time. This motion happens in phys- 
ical time t and the pdf for the operational time having value t* is (as density 
in t^, evolving in physical time t) given by (5.20). In fact, by substituting 

y = tt* we find 

/'OO /'OO 

/ gp(U,t)dU= g p (t,Qdt = l, Vt>0. (5.21) 

The operational time i* stands in analogy to the counting index n in Eqs. 
(3.5) and (5.4). The other process is the process described by Eq. (5.13), a 
spatial probability density for sojourn of the particle in point x evolving in 
operational time i*, 

up(x, U) = f a ,e(x, U) . (5.22) 

We get the solution to the Cauchy problem (2.1), namely the pdf u(x,t) = 
up(x,t) for sojourn in point x, evolving in physical time t, by averaging 
up(x,t^) with the weight function gp(t*,t) over the interval < t* < 00 
according to (5.19). 



-nXr' 1 



'sign k 



(5.16) 



U [s p + 



| K |a f 0signK 



(5.17) 



+ 



| K |a^signK 



(5.18) 
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6 Sample path for space-time fractional diffusion 



In the series representation (3.5) of the CTRW the running index n (the 
number of jumps having occurred up to physical time t) is a discrete op- 
erational time, proceeding in unit steps. To this index n corresponds the 
physical time t = t n , the sum of the first n waiting times, and in physical 
space the position x = x n , the sum of the first n jumps, see Section 3. 

Rescaling space and physical time by factor h and r, obeying the scaling 
relation 

^h a = \T p , (6.1) 
and introducing, by sending {h — > , r — > 0}, continuous operational time 

U~n\r p ~n^h a . (6.2) 

Then, in the series representation (3.5) we have two discrete Markov pro- 
cesses (discrete in operation time n), namely a random walk in the space 
variable x, with jumps X n , and another random walk (only in positive di- 
rection) of the physical time t, making a forward jump T n at every instant 
n. 

In the diffusion limit the spatial process becomes an a-stable process for 
the position x = x = x(t*), whereas the unilateral time process becomes 
a unilateral (positively directed) /3-stable process for the physical time t = 
i = i(t*). A sample path of a diffusing particle in physical coordinates can 
be produced by combining in the (t, x) plane the two random functions 

x x x(i*) , , . 

t = i=t(u), [ ' 

both evolving in operational time i*, both being Markovian and obeying 
stochastic differential equations 

dx = d(Levy noise of order a and skewness#) , , . 

dt = e?(one sided Levy noise of order j3) . " 

This gives us in the (t, x) plane the t* - parametrized particle path, and by 
elimination of t* we get it as x = x(t). We suggest to call this procedure 
"construction of a particle path by parametric subordination" . Note that 
the process t = T(t*) yielding the second random function in (6.3) has the 
properties of a subordinator in the sense of Definition 21.4 in |45j . 
Concerning notation : It is good to make a conceptual distinction between 
the position x of an individual particle and the variable x, likewise between 
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the physical time position i and the physical time variable t. When there 
are many particles we have overall densities for them and for these densities 
fractional diffusion equations. The pdf for the particle being in point x = x 
at operational time £*, that we denote by up(x, £*) = f a ,e(%, t*), satisfies the 
evolution equation (Eq. (5.13) re-written with un) 

d 

— up{x,U)= x D'$ui3(x,t*) y u(x,0)=5(x). (6.5) 

The pdf for the physical time being in t = t at operational time i*, that we 
denote by v(t, t*) = gp{t,t*), obeys the skewed fractional equation 

J-u(M*) = tD?.pv(t,U), v(t,0)=5(t). (6.6) 

Remind : In operational time two Markovian random functions x(t*), t(i*) 
occur, as random processes, individually for each particle. In physical coor- 
dinates we have the t*-parametrized random path described by (6.3). 
Remark : It is instructive to see what happens for the limiting value (3 = 1. 
In this case the Laplace transform of gp(t,t*) = g~\(t,t*) is exp(— t*s), im- 
plying gi(t, i*) = 5(t — t*), the delta density concentrated on i = t*. So, in 
this case, t = i*, operational time and physical time coincide. 



7 Numerical results 

In this Section, after describing the numerical schemes adopted, we shall 
show the sample paths for four case studies of symmetric (9 = 0) frac- 
tional diffusion processes: {a = 2, (3 = 0.90}, {a = 2, (3 = 0.80}, 
{q = 1.5, (3 = 0.90}, {q = 1.5, (3 = 0.80}. As explained in the previous Sec- 
tions, for each case we need to construct the sample paths for three distinct 
processes, the parent process x = Y(t*), the leading process t = T(t*) (both 
in the operational time) and, finally, the subordinated process x = X(t), 
corresponding to the required fractional diffusion process. For this purpose 
we proceed as follows for the required three steps. 

First, let the operational time t* assume N discrete equidistant values in 
a given interval [0, T], that is t* jn = nT/N, n = 0, 1, . . . , N. As a working 
choice we take T = 1 and ./V = 10 6 . Then produce iV independent identically 
distributed (iid) random deviates, Y%, Y2, ■ ■ ■ , Yn having a symmetric stable 
probability distribution of order a, see the book by Janicki [23] for a useful 
and efficient method to do that. Now, with the points 

n 

x = 0, x n = '^2x k , n=l,...,N , (7.1) 
k=l 
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the couples (t* t7l ,x n ), plotted in the (t*,x) plane (operational time, physical 
space) can be considered as points of a true sample path {x(t*) : < i* < T} 
of a symmetric Levy motion with order a corresponding to the integer values 
of operational time = In this identification of t* with n we use the 
fact that our stable laws for waiting times and jumps imply A = (i = 1 in 
the asymptotics (5.1) and (5.2) and r = h = 1 as initial scaling factors in 
(5.3) and (5.14). 

In order to complete the sample path we agree to connect every two suc- 
cessive points (t* :n ,x n ) and (t* >n+ i, x n+ \) by a horizontal line from (i* jn ,x n ) 
to (t*, n +i, x n ), and a vertical line from (t* in+ i, x n ) to (t*, ra +i, x n+ i). Ob- 
viously, this is not the 'true' Levy motion from point {t* )n ,x n ) to point 
(t* jn+ i, x n+ i), but from the theory of CTRW we know this kind of discrete 
random process to converge in the appropriate sense to Levy motion. The 
points (t*, n ,x n ) are points of a true Levy motion. 

As a second step, we produce N iid random deviates, T±, T2, . . . , T/v hav- 
ing a stable probability distribution with order and skewness — (3 (extremal 
stable distributions). Then, consider the points 

n 

to = 0, t n = Y,T k , n = l,...,N, (7.2) 
k=i 

and plot the couples (t* tn ,t n ) in the (t*,i) (operational time, physical time) 
plane. By connecting points with horizontal and vertical lines we get sample 
paths {£(£*) : < < Nt = 1} describing the evolution of the physical 
time t with increasing operational time t*. 

The final (third) step consists in plotting points (t{t if ^ n ),x{t 1f , n )) in the 
(t, x) plane, namely the physical time-space plane, and connecting them 
as before. So one gets a good approximation of the sample paths of the 
subordinated fractional diffusion process of parameters a, (3 and = 0. 

Now as the successive values of t* jn and x n are generated by successively 
adding the relevant standardized stable random deviates, the obtained sets 
of points in the three coordinate planes: (t*,t), (t*,x), (t,x) can, in view 
of infinite divisibility and self-similarity of the stable probability distribu- 
tions, be considered as snapshots of the corresponding true random processes 
occurring in continuous operational time and physical time t, correspond- 
ingly. Clearly, fine details between successive points are missing. They are 
hidden: 

- In the (i*, x) plane in the horizontal lines from (t*^, x n ) to (£*, n +:b x n ) and 
the vertical lines from (i^n+i, x n ) to (i^n+i, x n+ \). 

- In the (i*,t) plane in the horizontal lines from (t* jn ,t n ) to (i^n+i, t n ) and 
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the vertical lines from (i* jn +l>in) to (i* )n +i, tn+i). 

- In the (t, a;) plane in the horizontal lines from (t n ,x n ) to (t n +i,x n ) and 
the vertical lines from (t n+ i,x n ) to (t n+ i,x n+ i). 

The well-scaled passage to the diffusion limit here consists simply in reg- 
ularly subdividing the {i*} intervals of length 1 into smaller and smaller 
subintervals (all of equal length r and adjusting the random increments of t 
and x according to the requirement of self-similarity, namely taking, respec- 
tively, the waiting times and spatial jumps as r 1///3 multiplied by a standard 
extreme /3-stable deviate, t 1//q multiplied by a standard (in our special case: 
symmetric) o-stable deviate, respectively, as required by the self-similarity 
properties of the stable probability distributions) . Furthermore if we watch 
sample path in a large interval of operational time t*, the points (^* 1 n?^n) 
and (t*,n+l)^n+l) wu l m the graphs appear very near to each other in op- 
erational time t* and aside from missing mutually cancelling jumps up and 
down (extremely near to each other) we have a good picture of the true 
processes. 

The resulting sample paths for all the processes involved in the two case 
studies are presented in the Figs. 1-6. The figure captions should clarify our 
strategy. Figs. 1 and 4 are referring to the parent processes characterized 
by the parameter a = 2 and a = 1.5. Figs. 2 and 5 are devoted to the 
leading processes^ characterized by the parameter /3 = 0.9 and (3 = 0.8 in 
the Right and Left plates, respectively. As a consequence Figures 2 and 5 
are identical, because are referring to the same processes. Finally, Figs. 3 
and 6 are devoted to the subordinated processes resulting from the previous 
parent and leading processes. Specifically in Fig. 3 the Left and Right plates 
show sample paths for a = 2 and (3 = 0.9, 0.8, respectively, and in Fig. 6 
the Left and Right plates show sample paths for a = 1.5 and = 0.9,0.8, 
respectively. 

By observing the figures the reader will note that horizontal segments 
(waiting times) in the (t, x) plane (Fig. 3, Fig. 6) correspond to vertical 
segments (jumps) in the (t, i*) plane (Fig. 2, Fig. 5). Actually, the graphs in 
the (i, x)-plane depict continuous time random walks with waiting times Tj. 
(shown as horizontal segments) and jumps (shown as vertical segments). 
The left endpoints of the horizontal segments can be considered as snapshots 
of the true particle path (the true random process to be simulated), the 

1 Figs. 2 and 5 alternatively can also be viewed as graphical representations of the 
directing processes i, = T,(t) in the sense of Feller, see Section 4. We note that the 
directing processes, exhibiting horizontal segments, are no longer Levy processes even if 
the random functions i» = T*(t) are non-decreasing and right-continuous like the leading 
processes t = T(T*). This explains the non-Markovianity of the subordinated processes. 
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Figure 3: A sample path for the subordinated process x = X(t). 
LEFT: {a = 2, = 0.90}, RIGHT: {a = 2 , = 0.80}. 
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Figure 4: A sample path for the parent process x = Y(t*) with {a = 1.5}. 




Figure 5: A sample path for the leading process t = T(t*). 
LEFT: {p = 0.90}, RIGHT: {(3 = 0.80}. 




Figure 6: A sample path for the subordinated process x = X(t). 
LEFT: {a = 1.5, (3 = 0.90}, RIGHT: {a = 1.5, (3 = 0.80}. 
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segments being segments of our ignorance. In the interval t n < t < t n+ i the 
true process (namely the spatial variable x = X{t)) may jump up and down 
(infinitely) often, the sum (or integral) of all these ups (counted positive) 
and downs (counted negative) amounting to the vertical jump X n+ \. 

Finer details will become visible by choosing in the operational time t* 
the step length r smaller and smaller. In the graphs we can clearly see 
what happens for finer and finer discretization of the operational time t*, 
by adopting 10 1 , 10 2 , 10 3 steps, see Figures 7-18. As a matter of fact there 
is no visible difference in the transition for the successive decades 10 4 , 10 5 , 
10 6 steps as the great majority of spatial jumps and waiting times are very 
small. This property also explains the visible persistence of large jumps and 
waiting times even of very small steps r of the operational time. 

8 Conclusions 

Starting from the series representation (3.5) of the CTRW, by considering 
there the running index of summation as discrete operational time and pass- 
ing to the diffusion limit in a well-scaled way, we have shown how to arrive at 
the integral formula of subordination in fractional diffusion. Furthermore, 
we have explained how, in analogy to the construction of particle paths in 
CTRW, particle paths in space-time fractional diffusion can be obtained by 
composition of two stable (hence Markovian) processes (one for the physical 
time, the other for the position in space, both processes running in opera- 
tional time). By this composition we get in physical space-time the particle 
path, parametrized by the operational time. For this construction of a parti- 
cle path we suggest the name parametric subordination. The essential games 
are played in operational time, for construction of a particle path we avoid 
to explicitly run the hitting time process (see [32J ) generating from physical 
time the operational time. 
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Figure 9: A sample path for the parent process x = Y(t*). 
LEFT: {a = 2, N = 10 3 }, RIGHT: {a = 1.5, N = 10 3 }. 
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Figure 10: A sample path for the leading process t = T(t*). 
LEFT: {0 = 0.9, N = 10 1 }, RIGHT: {(3 = 0.8, N = 10 1 }. 
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Figure 11: A sample path for the leading process t = T(t*). 
LEFT: {(5 = 0.9, N = 10 2 }, RIGHT: {(3 = 0.8, N = 10 2 }. 
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Figure 12: A sample path for the leading process t = T(t*). 
LEFT: {/3 = 0.9, N = 10 3 }, RIGHT: {(3 = 0.8, N = 10 3 }. 
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Figure 13: A sample path for the subordinated process x = X(t). 
LEFT: {a = 2, (3 = 0.90, TV = 10 1 }, RIGHT: {a = 2, /3 = 0.80, TV = 10 1 }. 
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Figure 14: A sample path for the subordinated process x = X(t). 
LEFT: {a = 2, /3 = 0.90, TV = 10 2 }, RIGHT: {a = 2, (3 = 0.80, TV = 10 2 }. 





Figure 15: A sample path for the subordinated process x = X(t). 
LEFT: {a = 2, (3 = 0.90, /V = 10 3 }, RIGHT: {a = 2, /? = 0.80, TV = 10 3 }. 
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Figure 16: A sample path for the subordinated process x = X(t). 
LEFT: {a = 1.5, = 0.90, N = 10 1 }, RIGHT: {a = 1.5, = 0.80, N = 10 1 }. 




t t 



Figure 17: A sample path for the subordinated process x = X(t). 
LEFT: {a = 1.5, = 0.90, N = 10 2 }, RIGHT: {a = 1.5, = 0.80, A" = 10 2 }. 





Figure 18: A sample path for the subordinated process x = X(t). 
LEFT: {a = 1.5, = 0.90, N = 10 3 }, RIGHT: {a = 1.5, = 0.80, N = 10 3 }. 
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